#include "newton.h"

#include <cmath>

namespace NumericalAnalysis
{
    // Performs the Newton method to estimate a value.
    // It takes a NewtonInputData and returns a NewtonOutputData.
    NewtonOutputData NewtonMethod(const NewtonInputData& inputData)
    {
    	assert(inputData.mEstimationFunction && "NewtonMethod: NULL pointer");
    	assert(inputData.mDerived && "NewtonMethod: NULL pointer");

    	// Instantiate the output data.
    	NewtonOutputData outputData;
    	
    	// While we must continue the iterations 
    	// and we do not reach the max number of iterations,
    	// we will perform another iteration.
    	double previousEstimation = 0.0;
    	double currentEstimation = inputData.mInitialEstimation;
    	uint32_t currentIteration = 0;
    	bool continueIterations = true;
    	while(continueIterations && currentIteration < inputData.mMaxIterations)
    	{
    		// Increment the current iteration.
    		++currentIteration;

    		// Check if the absolute value of the current estimation 
    		// valued by the function is enoughly near
    		// using the stop bound or the current relative error
    		// is enoughly small according to the stop bound or 
    		// if the derive of the estimation function, valued on the current estimation
    		// is zero. If this is the case, we will stop the iterations.
    		const double valuedCurrentEstimation = inputData.mEstimationFunction(currentEstimation); // f(Xn)
    		
    		const double result = (currentEstimation - previousEstimation) / currentEstimation; // (Xn - Xn-1) / Xn
    		const double relativeError = abs(result); // | (Xn - Xn-1) / Xn |
    		const double derivedCurrentEstimation = inputData.mDerived(currentEstimation); // f'(Xn)

            const double absValuedCurrentEstimation = abs(valuedCurrentEstimation); // | f(Xn) |
    		if(absValuedCurrentEstimation < inputData.mStopBound 
				|| relativeError < inputData.mStopBound 
				|| derivedCurrentEstimation == 0.0)
    		{
    			// Stop the iterations.
    			continueIterations = false;

    			// Set the final iteration.
    			outputData.mFinalIterationNumber = currentIteration;

    			// Set the final approximation.
    			outputData.mFinalEstimation = currentEstimation;

    			// Set the absolute value of the final approximation 
    			// valued by the function.
    			outputData.mAbsOfValuedFinalEstimationByFunction = absValuedCurrentEstimation;

    			// Set the relative error.
    			outputData.mRelativeError = relativeError;
    		}

    		// Else update the current estimation.
    		else
    		{
    			// Calculate the next estimation.
    			const double nextEstimation = currentEstimation - (valuedCurrentEstimation / derivedCurrentEstimation); // Xn - (f(Xn) / f'(Xn))

    			// Update the previous estimation.
    			previousEstimation = currentEstimation;

    			// Update the current estimation.
    			currentEstimation = nextEstimation;
    		}
    	}

    	// Return the output data.
    	return outputData;
    }
}
